Ejemplo de análisis espacial: densidad de la red vial

Ejemplo de análisis espacial: densidad de la red vial#

Abrir en Colab

Introducción#

La densidad de la red vial para un polígono se define como:

\[longitudRedVial(km) / áreaPolígono(km^{2})\]

Por ejemplo, si un cantón tiene 500 km de longitud de red vial y un área de 1000 km2, la densidad de su red vial es 0.5.

Este cuaderno de notas calcula la densidad de la red vial para cada uno de los cantones de Costa Rica y la presenta en un mapa de coropletas.

Carga de bibliotecas#

# Carga de geopandas con el alias gdp
import geopandas as gpd

# Otras bibliotecas

# Carga de pandas con el alias pd
import pandas as pd

# Carga del módulo pyplot de matplotlib con el alias plt
import matplotlib.pyplot as plt

# Carga de la clase WebFeatureService del módulo wfs de owslib
# Permite interactuar con servicios web geoespaciales tipo WFS
from owslib.wfs import WebFeatureService

# Carga de la clase BytesIO del módulo estándar io
# Permite crear un objeto en memoria que actúa como un archivo binario
from io import BytesIO

# Carga de la biblioteca Folium, para mapas interactivos
import folium

Carga de datos#

Cantones#

# Conexión al servicio WFS del IGN 5k CO
wfs_url = 'https://geos.snitcr.go.cr/be/IGN_5_CO/wfs'
wfs_version = '1.1.0'
wfs = WebFeatureService(url=wfs_url, version=wfs_version)

# Obtener la capa de cantones
capa = 'IGN_5_CO:limitecantonal_5k'
respuesta = wfs.getfeature(typename=capa, outputFormat='application/json')

# Leer la respuesta en un GeoDataFrame
cantones_gdf = gpd.read_file(BytesIO(respuesta.read()))

# Reducción de columnas
cantones_gdf =cantones_gdf[['CÓDIGO_CANTÓN', 'CANTÓN', 'geometry']]

# Definir un índice
# cantones_gdf.set_index('CÓDIGO_CANTÓN', inplace=True)

# Mostrar las primeras filas del GeoDataFrame
print(cantones_gdf.head())

# Mapa
cantones_gdf.plot()

   CÓDIGO_CANTÓN        CANTÓN  \
0            101      San José   
1            102        Escazú   
2            103  Desamparados   
3            104      Puriscal   
4            105       Tarrazú   

                                            geometry  
0  MULTIPOLYGON (((481030.328 1102633.167, 481032...  
1  MULTIPOLYGON (((479844.919 1102669.584, 479851...  
2  MULTIPOLYGON (((494091.231 1094936.82, 494091....  
3  MULTIPOLYGON (((455934.039 1095770.324, 455949...  
4  MULTIPOLYGON (((501998.931 1074558.189, 502000...  

<Axes: >
../../_images/8ff5b6d5e88c6b868a6e9000873c37782fc09c6e667f849edfccc463632384ee.png

Red vial#

# Conexión al servicio WFS del IGN 200 k
wfs_url = 'https://geos.snitcr.go.cr/be/IGN_200/wfs'
wfs_version = '1.1.0'
wfs = WebFeatureService(url=wfs_url, version=wfs_version)

# Obtener la capa de red vial
capa = 'IGN_200:redvial_200k'
respuesta = wfs.getfeature(typename=capa, outputFormat='application/json')

# Leer la respuesta en un GeoDataFrame
red_vial_gdf = gpd.read_file(BytesIO(respuesta.read()))

# Añadir una columna con un identificador único
red_vial_gdf['red_vial_id'] = range(1, len(red_vial_gdf) + 1)

# Reducción de columnas
red_vial_gdf =red_vial_gdf[['red_vial_id', 'geometry']]

# Definir un índice
# red_vial_gdf.set_index('red_vial_id', inplace=True)

# Mostrar las primeras filas del GeoDataFrame
print(red_vial_gdf.head())

# Mapa
red_vial_gdf.plot()

/home/mfvargas/miniconda3/envs/gf0657-2024-ii/lib/python3.12/site-packages/pyogrio/raw.py:196: RuntimeWarning: Several features with id = 0 have been found. Altering it to be unique. This warning will not be emitted anymore for this layer
  return ogr_read(

   red_vial_id                                           geometry
0            1  LINESTRING (543973.13 1157346.96, 543972.52 11...
1            2  LINESTRING (542005.8 1160843.73, 542304.68 116...
2            3  LINESTRING (542872.85 1160705.42, 542946.13 11...
3            4  LINESTRING (527814.53 1167099.79, 527692.44 11...
4            5  LINESTRING (538428.65 1157738.05, 538251.35 11...

<Axes: >
../../_images/2111dbcdc954484387a40a7f293848547d5adc021d61876398cdc2901c89acab.png

Análisis#

# Generar un segmento de red vial 
# por cada intersección con un polígono de cantón
interseccion_redvial_cantones_gdf = gpd.overlay(
    red_vial_gdf, 
    cantones_gdf, 
    how='intersection'
)

interseccion_redvial_cantones_gdf.head()

/tmp/ipykernel_228903/3535146563.py:3: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: EPSG:5367
Right CRS: EPSG:8908

  interseccion_redvial_cantones_gdf = gpd.overlay(
red_vial_id CÓDIGO_CANTÓN CANTÓN geometry
0 1 702 Pococí LINESTRING (543973.13 1157346.96, 543972.52 11...
1 2 702 Pococí LINESTRING (542005.8 1160843.73, 542304.68 116...
2 3 702 Pococí LINESTRING (542872.85 1160705.42, 542946.13 11...
3 4 702 Pococí LINESTRING (527814.53 1167099.79, 527692.44 11...
4 5 702 Pococí LINESTRING (538428.65 1157738.05, 538251.35 11... 
# Calcular la longitud de cada segmento de red vial intersecado
interseccion_redvial_cantones_gdf["LONGITUD_RED_VIAL"] = interseccion_redvial_cantones_gdf["geometry"].length/1000

interseccion_redvial_cantones_gdf.head()
red_vial_id CÓDIGO_CANTÓN CANTÓN geometry LONGITUD_RED_VIAL
0 1 702 Pococí LINESTRING (543973.13 1157346.96, 543972.52 11... 8.744950
1 2 702 Pococí LINESTRING (542005.8 1160843.73, 542304.68 116... 8.188802
2 3 702 Pococí LINESTRING (542872.85 1160705.42, 542946.13 11... 2.135535
3 4 702 Pococí LINESTRING (527814.53 1167099.79, 527692.44 11... 2.302054
4 5 702 Pococí LINESTRING (538428.65 1157738.05, 538251.35 11... 8.779122 
# Agrupación y suma de los los segmentos de red vial por cantón
suma_redvial_cantones_df = interseccion_redvial_cantones_gdf.groupby(["CÓDIGO_CANTÓN"])["LONGITUD_RED_VIAL"].sum()

# Convertir la serie en dataframe
suma_redvial_cantones_df = suma_redvial_cantones_df.reset_index()

# Desplegar longitudes de red vial por cantón
suma_redvial_cantones_df.sort_values(by='LONGITUD_RED_VIAL', ascending=False)
CÓDIGO_CANTÓN LONGITUD_RED_VIAL
29 210 2591.051204
18 119 1429.548737
65 601 1319.628280
67 603 1297.054690
53 410 1103.343280
... ... ...
14 115 23.190868
50 407 20.224896
12 113 16.403583
52 409 15.086709
51 408 11.006955

84 rows × 2 columns

# Agregar a la capa de cantones la longitud de la red vial para cada polígono
cantones_redvial_gdf = cantones_gdf.merge(suma_redvial_cantones_df, on="CÓDIGO_CANTÓN")

# Calcular el área del cantón
cantones_redvial_gdf["ÁREA_CANTÓN"] = cantones_redvial_gdf["geometry"].area/1000000

# Calcular la densidad de la red vial
cantones_redvial_gdf["DENSIDAD_RED_VIAL"] = cantones_redvial_gdf["LONGITUD_RED_VIAL"] / cantones_redvial_gdf["ÁREA_CANTÓN"]

# Desplegar la densidad de red vial en cada cantón
cantones_redvial_gdf[['CÓDIGO_CANTÓN', 'CANTÓN', 'DENSIDAD_RED_VIAL', 'ÁREA_CANTÓN', 'LONGITUD_RED_VIAL']].sort_values(by="DENSIDAD_RED_VIAL", ascending=False)
CÓDIGO_CANTÓN CANTÓN DENSIDAD_RED_VIAL ÁREA_CANTÓN LONGITUD_RED_VIAL
0 101 San José 2.259410 44.624284 100.824537
7 108 Goicoechea 2.148213 31.696493 68.090826
17 118 Curridabat 2.065924 16.062606 33.184129
12 113 Tibás 1.983644 8.269421 16.403583
19 120 León Cortés Castro 1.954404 121.894583 238.231301
... ... ... ... ... ...
69 605 Osa 0.351286 1932.697429 678.930001
78 701 Limón 0.291557 1769.378451 515.874969
77 613 Puerto Jiménez 0.287517 720.432185 207.136511
44 401 Heredia 0.287120 283.109137 81.286423
81 704 Talamanca 0.117583 2792.225947 328.318153

84 rows × 5 columns

# Mapa estático
cantones_redvial_gdf.plot(
    column="DENSIDAD_RED_VIAL", 
    legend=True,
    cmap='OrRd', 
    scheme='quantiles',
    figsize=(20, 20)
)

<Axes: >
../../_images/6172aeb899652ec2e78bc47db3faff0a42025c0cd7ad3f25e007cfedc8354d6e.png 
# Antes de agregar cantones_redvial_gdf, 
# se simplifica (i.e. se reduce el número de vértices)
# para generar un HTML más pequeño (< 100 MB)

# Hacer una copia de cantones_redvial_gdf para simplificarla
cantones_redvial_simplificado_gdf = cantones_redvial_gdf

# Definir el factor de simplificación
# Cuanto mayor sea el valor, mayor será la simplificación
factor_simplificacion = 0.01

# Simplificar geometría
cantones_redvial_simplificado_gdf['geometry'] = cantones_redvial_simplificado_gdf['geometry'].simplify(
    factor_simplificacion, 
    preserve_topology=True
)

# Crear el mapa interactivo con la densidad de la red vial
m = cantones_redvial_simplificado_gdf.explore(
    column='DENSIDAD_RED_VIAL',
    name='Densidad de la red vial en cantones',
    cmap='OrRd',
    tooltip=['CANTÓN', 'DENSIDAD_RED_VIAL'],
    legend=True,
    legend_kwds={
        'caption': "Densidad de la red vial", 
        'orientation': "horizontal"
    },
)

# Añadir la red vial al mapa
#red_vial_gdf.explore(
#    m=m, # se usa el mapa que se creó en la instrucción anterior
#    name='Red vial',
#    popup=True
#)

# Agregar un control de capas al mapa
folium.LayerControl().add_to(m)

# Mostrar el mapa interactivo
m
Make this Notebook Trusted to load map: File -> Trust Notebook

Ejercicios#

  1. Programe un mapa interactivo de coropletas que muestre la densidad de la red vial en las regiones socioeconómicas de Costa Rica.

Contents